Method for automatic measurement of the doppler angle and arrangement for carrying out the method

ABSTRACT

Automatic measurement is made of the angle DA enclosed by the direction of ultrasonic echographic excitation and the axis of a vessel in an echographic image on the basis of prior designation of an initial point P i  in the vessel. A first isotropic tracing of rays starting from the initial point P i  is used to produce a histogram of the grey levels of the points of the rays. An algorithm is then applied to the histogram in order to classify the grey levels of selected points. A second tracing of rays is made which is restricted to the walls of the vessel and results in a local mark with triangular sectors from which the slope (a) of a regression line is determined and the calculation of the Doppler angle DA is made as DA=Arc tg(a).

The invention relates to a method for automatic measurement of theDoppler angle DA, enclosed by the direction of an echographic excitationand the axis of a blood vessel in an ultrasonic echographic grey-levelimage, on the basis of prior designation of an initial point situated inthe vicinity of the axis of the relevant vessel.

The invention also relates to an arrangement for carrying out themethod.

In order to provide a quantitative measure of the speed of blood inarteries, a method for the ultrasonic imaging of speeds by colourencoding is based on knowing the angle between the incident ultrasonicbeam and the local axis of the vessel. This holds for all ultrasonicsignal processing methods currently carried out to produce acolour-encoded representation of flows in any commercially availableechograph. Moreover, reliable determination of this angle, referred toas the Doppler angle, is necessary to extract echographic data ofparameters whose measurement is rather complex and to make furtherprogress in the study of the behaviour of arteries.

Considering the foregoing, it would be justified to conceive and carryout an automatic process for the extraction of the value of the Dopplerangle from a medical ultrasonic image with a high precision, for exampleof the order of 1 degree.

The Doppler angle still is a difficult aspect of all commerciallyavailable apparatus: the radiologist must determine the value thereof ina rather approximative manner by judging the alignment of a small linesegment traced in the image with respect to the axis of the blood vesselanalyzed.

From Japanese Patent No. 5-31110 it is known to detect automatically theangle enclosed by the direction of the blood flow with respect to thedirection of an echographic beam. Via a transmission using anechographic probe comprising several transducers, followed by a separateanalysis in respect of phase and frequency of the reflected beamsreceived in two receiving sections of the probe, a given Dopplerdisplacement is determined wherefrom a speed vector can be deduced atthe point affected by the probe; this corresponds to the determinationof a Doppler angle. However, it is very difficult to carry, out a directmeasurement without an intermediate echographic image, because it is notensured that the axis of the blood vessel probed is indeed present inthe plane of the image of the artery as is necessary for the measurementof the Doppler angle. When the axes of the echographic excitation andthat of the vessel are not situated in the same plane, the detection ofthe Doppler angle remains possible, but it is liable to be affected bysubstantial errors; this is not acceptable when the required precisionis of the order of approximately 1 degree.

It is an object of the invention to provide a method for automaticmeasurement of the Doppler angle of a trace of a blood vessel in anechographic image, starting from a designated point in the directvicinity of said trace.

It is another object of the invention to provide a method for automaticmeasurement of the Doppler angle which is exact and fast.

It is another object to provide an arrangement for carrying out theabove measuring method.

These objects are achieved, and the drawbacks of the prior art aremitigated, in that the method of the kind set forth is characterized inthat it includes the following steps:

a) a first isotropic tracing of rays from said initial point, said raystraversing the echographic image partly or entirely so as to provide ahistogram of the grey levels of selected points which are regularlydistributed along said rays,

b) executing an image processing algorithm which is applied to saidhistogram in order to classify the grey levels of said selected pointsinto at least two categories CL₁, CL₂, . . . , two adjacent classesbeing separated by thresholds T₁, T₂, . . . , expressed in grey levels,one of the classes CL_(i), bounded by the lower threshold T_(i-1), beingrepresentative of walls of blood vessels in said image,

c) a second tracing of rays from said initial point, during which thegrey level of each point of each ray is compared with the thresholdT_(i-1) and each ray is restricted to the first end point Pl encounteredwhose grey level is equal to or higher than T_(i-1), resulting in arepresentation of said blood vessel in the form of a first local markwhich is composed of triangular sectors which have said starting pointwithin said first local mark in common,

d) determination of the slope a of the regression line of said firstlocal mark by application of a linear regression method to the N pixels,having the coordinates x(n) and y(n) of said first local mark, as:##EQU1## where: ##EQU2## e) and calculation of the Doppler angle as:

    DA=Arc Tg (a).

As is often the case, an echographic image may contain the trace ofseveral vessels. Therefore, the radiologist wishing to measure the angleshould first designate the vessel of interest, for example by way of themouse of a workstation on the screen of which the echographic imagewould be displayed. Said initial point can thus be designated, afterwhich the method for measuring the Doppler angle in accordance with theinvention is automatically performed in an exact and very fast manner.

An attractive version of the invention is characterized in that themethod described above includes the supplementary step of validationwhich consists in validating the value of the Doppler angle DA only fora correlation coefficient r of the pixels of said first local mark whichexceeds a threshold Rm of predetermined value between 0 and 1, where:##EQU3## The value chosen for Rm is preferably larger than 0.5, thusimposing a minimum quality criterion on the angle measurement performed.

For the calculation of the slope a in conformity with the expression (1)and of the coefficient r in conformity with the expression (2) in afurther preferred version of the invention it is particularly attractiveto perform the calculation of the terms involved in the expressions (1)and (2) for each triangular sector, formed by two rays whose ends havethe coordinates x(m), y(xm) and x(M), y(xM) (where XM≧xm), on the basisof the following analytical formula: ##EQU4## in which the functionf(x,y) is given the following successive values: 1; x; y; xy; xx; yy,and the terminations x₀, x₁, y₀, y₁, x₂, y₂, x₃ are simple linearfunctions of coordinates of said end points of the two rays (notably: x₀=0, x₁ =xm, x₂ =xM), with 6 sets of possible values associated with 6respective classes of different triangular sectors.

It is to be noted that in the above method for automatic measurement inthe step d) for determining the slope a, being representative of theslope of the axis of said blood vessel, this slope a can be obtained as:##EQU5## in which a is the slope of the inertial axis of said firstlocal mark.

Another version can also be envisaged for the supplementary step forvalidating the value of the Doppler angle DA automatically measured bymeans of the method in accordance with the invention. This supplementaryvalidation step consists, for example, in validating the Doppler angleDA only if the symmetry of said local mark with respect to itsregression line (or its inertial axis) is considered to be sufficient,said symmetry being deduced from the calculation of the centred momentof order 3, being sk(f): ##EQU6## and f is the function:

    y(n)=f(x(n)), where n is within  1, . . . , N!,

the centred moment of order 3 being calculated for two curves which arerepresentative of the function f:

sk1 for the curve formed by the extremities of rays of said local markwhich are situated above the regression line (or the inertial axis)taken as the axis of the abscissae (y>0),

sk2 for the curve formed by the extremities of rays situated below saidline (y<0),

the values of sk1 and sk2 being compared afterwards.

In order to perform the comparison, a start is made from the mark of theaxis determined by the slope a, previously defined, and first sk1 andsk2 are calculated as indicated above. Subsequently, the following testsare successively carried out in order to obtain the reliabilitycoefficient:

1--if the absolute value of sk1 and absolute value of sk2 are smallerthan 1, the mark is reasonably symmetrical:

    reliability=100%

2--if not, if sk1 and sk2 have the same sign (comparable dissymmetry ofthe two curves), the following estimate is made:

    reliability=100.0*(1-|sk1-sk2|/|sk1+sk2.vertline.)%

3--if not, the mark is not sufficiently symmetrical with respect to itsaxis:

    reliability=0%

In order to ensure that the automatic measurement of the Doppler angleis correct, the radiologist should select the initial point within thevessel for which the Doppler angle is to be measured and, even if it issituated therein, it may occur that the initial point chosen is situatedvery close to the wall of the vessel; this is not an optimum situation.Another preferred version of the invention, enabling further enhancementof the measurement of the Doppler angle, is characterized in that itincludes a supplementary step for positioning which consists incalculating the centre of gravity G of said first local mark, saidcentre of gravity G then being chosen as the new initial point which hasbeen optimized for optimized determination of the Doppler angle DA byway of the method for measurement as defined in the foregoingparagraphs. This operation, consisting in shifting the initial point toa new initial point which is formed by the centre of gravity of thepreviously obtained local mark, can advantageously be repeated severaltimes.

In order to carry out the method of the invention, it is proposed to usean arrangement (preferably a special-purpose workstation) whichincludes:

a memory for storing grey-level values of the echographic image in theform of a matrix of pixel values,

a screen displaying the echographic image in the form of atwo-dimensional matrix of pixels,

pointer means for said initial point in the image, and

calculation means having access to said memory and to the pixel which isrepresentative of said initial point, said calculation means beingspecifically programmed so as to execute the algorithms and calculationsnecessary for carrying out the method.

The automatic and exact supply of the Doppler angle as permitted by thepresent invention facilitates the interpretation of given information ofultrasonic origin on the screen of the echography apparatus, or on thatof a workstation connected thereto, and makes it more interactive forthe operator.

All ultrasonic echography systems available at present utilize themeasurement of the speed of a blood flow via one and the same techniquewhich is based on the Doppler effect. The differences between systemsoccur at the level of the implantation for which two techniques coexist:the so-called "frequential" approach of the conventional Doppler and the"temporal" approach developed notably by Philips; in this respectreference is made to European Patent No. 0 225 667 (PHF 85/593).

In all practical cases ultrasonic pulses are injected into the medium tobe analyzed, presumably containing a blood vessel, with a period T(Pulse Repetition Frequency or PRF) and the ultrasonic signal returnedby the medium in response to each pulse is recorded by the echograph(lines RF). Each sample of the lines RF thus represents the reflectivityof the medium at a given depth z or, equivalently, corresponds to agiven time of flight t for the pulse. It is to be noted that the signalreturned has travelled from the probe to the scatter centre and back; ifc is the speed of sound in the medium analyzed, therefore:

    z=ct/2.

In a window situated at a given depth on a line RF (for example, in ablood vessel), the signal slice corresponds to a "signature" whichbelongs to the instantaneous spatial distribution of the scatter centres(blood cells) in the interacton volume (resolution cell) between themedium analyzed and the ultrasonic beam. This signature belongs to aparticular disposition of scatter centres in the resolution cell,because it is the result of interference between all secondary sourcesconstituting these scatter centres (speckle); statistically speaking itis quasi inimitable.

Consequently, if the PRF is low enough so that the scatter centres donot have time to change their spatial disposition significantly betweentwo successive RF line acquisitions, a same group of scatter centreswill impress strongly resembling signatures from one line to another,and its shift in the depth direction can be detected and followed in thecourse of time. It is during this phase of following at a given depththat the temporal technique (referred to as CVI or Colour VelocityImaging) and the frequential (or Doppler) technique differ: the Dopplertechnique marks the shift of a signature while analyzing the evolutionof the local frequency of the successive lines RF, whereas the CVItechnique directly analyses the temporal signal of the lines RF.

The shift of a given signature between two successive lines enablesrecourse to the speed v of the group of associated scatter centres viathe time shift dt estimated by the CVI (dt=2vT/c) or the dephasingestimated by the Doppler technique (df=2πf dt, where f is the ultrasonicfrequency of the pulse). Thus, these two methods are equivalent intheory.

If the angle DA enclosed by the axis of the ultrasonic beam and thedirection of the flow of scatter centres is called the Doppler angle,the speed of the group of scatter centres considered, observed on thepropagation axis of the pulses, will be the projection v cos(DA) of itsreal speed v, it being readily possible to derive therefrom theexpression:

    v=c dt/(2T cos (DA)).

Consequently, the Doppler angle constitutes a fundamental parameter forthe estimation of the speed of a blood flow in an ultrasonic echograph,regardless of the method used for this estimation.

At present there are four major types of applications (referred toherein as A, B, C, D) enabling the encoding or usage of the speedinformation of a blood flow in a manner suitable for use in the medicalfield:

    ______________________________________                                        A      the representation of the speed information in the form                       of an image or colour speed map,                                       B      the representation in the course of the evolution time of                     the dynamic of the speeds or a histogram of the speed                         distribution (Doppler spectrum),                                       C      the determination of the blood flow rate (Philips CVI-Q),              D      the extraction of arterial parameters linked to the study of the              dilatation of the vessels in the course of the cardiac                 ______________________________________                                               cycle.                                                             

The first three of these techniques appear to be functions linked to themeasurement of the Doppler angle by way of the dependence of the speedon this parameter. The basic expressions for estimating the speed of agroup of scatter centres on the basis of a time shift dt (CVI) or afrequency shift df (Doppler) confirm:

CVI: v=cdt/(2T cos (DA))

Doppler: v=cdf/(4π fr cos (DA))

The technique D for that matter utilizes a dependency on sin(DA).

Principle of Colour Speed Mapping (Application A).

An echograph display, by default, a reflectivity image in grey levelswhich is formed by a juxtaposition of lines RF whose envelope has beentaken and which are recorded side-by-side in order to cover the zone ofinterest (mode B). These lines do not relate to the lines used forestimating the speed. A map of the speed distribution is formed bysuperposing on said reflectivity image a dynamic image (in the course ofthe cardiac cycle) whose colour scale corresponds to the scale of theinstantaneous speeds observed in conformity with a given correspondencerule (colour intensity proportional to the speed, colour a function ofthe sign of the speed, . . .).

In practice supplementary pulses (referred to as "colours") are emittedby the probe during the conventional image scanning and the aboveprocessing operation is applied thereto in order to detect shifts in theimage region of interest (presumably containing a blood vessel). Theparameter dt (or df) is estimated in all locations of the lines RF"colours" corresponding to a pixel of the region of interest, and acolour is assigned to said pixel in the dynamic image simply byapplication of the speed-colour correspondence rule (or equivalentlydt-colour or df-colour).

In these circumstances a version of the method of the invention, adaptedto CVI (temporal or frequential), comprises a supplementary step whichconsists of automatically correcting said speed colour scale so as toensure that it contains the dynamic of the speeds observed along theaxis of said blood vessel, by automatically displaying the maximumpositive and negative values of these speeds at the two extremities ofsaid scale.

Principle of the Speed Distribution Histogram (Application B)

An echographic image (grey levels or colour) corresponds to the emissionof pulses which are spatially shifted relative to one another in orderto scan the zone of interest of the medium. An alternative consists ofinjecting the successive pulses in the same location of the medium, sothat the evolution of the lines RF in the course of time (mode M) can beobserved. If "colour" excitations are mixed in this emission and areprocessed in the same way as described above or the estimation of theparameter dt (or df) at the locations of the blood flow, the evolutionof the speed distribution in the course of time can be observed in therelevant segment of the vessel.

The representation of this evolution takes the form of a histogramshowing the dynamic of the speeds as a function of time (either dt or dfas a function of time). The histogram in question, corresponding to atrue spectrum when it stems from the frequency analysis, is called aDoppler spectrum. The principle of its composition is very simple: theextremities of the dynamic of the parameter observed on thecorresponding colour lines are considered at the instant t before theyare represented, and there is obtained dt min<dt<dt max. Therefrom, thevalue of the Doppler angle DA being known, there is deduced the speeddynamic to be represented for this instant t:

    c dt min/(2T cos (DA))<v<c dt max/(2T cos (DA)).

Thus, this speed dynamic is represented by plotting a line vertically tothe abscissa t of the histogram, limited and governed by the extremespeed values observed.

In these circumstances an attractive version of the method for automaticmeasurement in accordance with the invention, in which the evolution ofthe speed distribution in the course of time in the relevant segment ofthe vessel in the mode M is represented, on a screen provided for thispurpose, in the form of a histogram showing the dynamic of the speeds asa function of time, referred to as a Doppler spectrum, comprises asupplementary step which consists of automatically correcting the scaleof the speeds represented in said Doppler spectrum in order to ensurethat it contains the dynamic of the speeds observed along the axis ofsaid blood vessel.

It is to be noted that for all envisaged applications of the inventionit is advantageous when the automatically measured value of the Dopplerangle DA is automatically displayed in a field of the screen containingthe echographic data.

These and other aspects of the invention will be apparent from andelucidated with reference to the embodiments described hereinafter.

In the drawings:

FIG. 1 shows the diagram of a workstation specifically arranged to carryout the invention.

FIG. 2 shows an echographic image in which an initial point is selected.

FIG. 3 illustrates the first tracing of rays from the initial point.

FIG. 4 shows the second tracing of rays, producing the first local mark.

FIG. 5 shows the regression line deduced from the first local mark.

FIG. 6 shows an organigram summarizing the execution of the entiremethod of the invention.

FIG. 7 shows the various feasible shapes of triangular sectors of thelocal mark.

For the sake of clarity, the contrast has been reversed in the FIGS. 2,3, 4, 5.

The workstation shown in FIG. 1 comprises, in known manner, amicroprocessor 1 which is connected to a memory 4, via a bidirectionaldata bus 2 and an address bus 3, and comprises two peripheral apparatuswhich are formed by a monitor 5 for the display of data stored in thememory 4 on a screen 6 and a mouse 7 which enables the marking andtracing of points on the screen 6. When an echographic image 8 whosedata are stored in the memory 4 in the form of pixels representative ofgrey levels is displayed on the screen 6, a radiologist is faced by thecustomary problem of determining the orientation of traces such as 9which stem from given blood vessels in the image whose echographicexcitation direction is represented by the vector 10. Therefore,normally a small segment of a straight line which is orientedapproximately in conformity with the estimated largest diameter of therelevant trace is superposed on the relevant trace by means of themouse. The angle searched, being the so-called Doppler angle which canbe readily calculated by means of the processor 1, is then the angleenclosed by said small line segment of line (not shown) with respect tothe vertical (in FIG. 1) which is identical to the direction of theechographic excitation (vector 10). This method is time-consuming,intricate and not very accurate. In order to improve this situation, themicroprocessor 1 is arranged to calculate said Doppler angle in a fastand exact manner on the basis of the designation of a single point,referred to as the initial point, in a trace of the blood vessel chosenby the user (the radiologist, the operator). To this end, aspecial-purpose processor 11, comprising the processor 1, also includescalculation means 12 necessary for carrying out the invention asdescribed hereinafter in several sections which successively describethe principle of the method and its execution details, and in an annexwhich specifies given calculations in detail.

Generally speaking, the method comprises two steps: segmentation of thelocal trace 9 of the vessel in the image 8 and linear regression overthe assembly of resultant points, referred to as a local mark, in orderto establish the direction of its axis.

An echographic image may contain the traces of several vessels which maybe of interest to a radiologist, so that it is necessary for theoperator to indicate the vessel of interest in the image before anytreatment can intervene, by means of the mouse of the system as shown inFIG. 2 in which the initial point Pi is selected within the local trace9 of a blood vessel.

Subsequently, in order to characterize alignment with respect to theshape of the vessel chosen, first its local mark is extracted from thesurrounding grey level in the image. This segmentation of the vessel isactually the detection of the walls of the vessel around the initialpoint indicated by the user (the practician, the physician).

Therefore, the user must perform the pointing operation in a uniformdark part of the vessel and it is also necessary to have a thresholdvalue available (referred to as T_(i-1)) which separates thereflectivity levels of the wall 11bis and the blood reflectivity levels(trace 9), thus leaving the operator a given margin so as to avoid smalltissue pans which are hardly visible and which could be present insidethe vessel.

For the determination of such a threshold value, preferably made of aconventional algorithm which is called ISODATA in the work: "Traitementnumerique des images", by MURAT KUNT, Collection Technique etScientifique des Telecommunications Presses Polytechniques etUniversitaires Romanes. This algorithm is based on a histogram of greylevels of pixels which comprises at least two classes. Preferably, ahistogram comprising three classes is chosen: CL₁ =VESSEL (lowreflectivity), CL₂ =INTERMEDIATE and CL₃ =WALL (high reflectivity). Thethreshold value T₂ which lies between the class INTERMEDIATE and theclass WALL is used in order to achieve the desired reliable walldetection.

The thresholds enable determination as to whether the starting pointsuitably belongs to the class VESSEL; this may urge the operator tostart pointing anew if the initial location is erroneous, thusconstituting an a posteriori validation of the initial point.

For as fast as possible operation, use is made of a partial histogramwhich utilizes a ray tracing method: as from the chosen initial pointPi, a predetermined number of radial rays which are angularlydistributed in an isotropic fashion, are traced in the image until theedge is reached, and only the pixels encountered along rays participatein the formation of the histogram. This ray tracing method is shown inFIG. 3.

When the ISODATA algorithm has produced the required threshold value(T_(i-1)), the same ray tracing method as described above is used tofind the first pixel P1 relating to a wall along each ray (it is ensuredthat a start is made within the vessel because of prior intervention bythe user): progression along each ray stops as soon as the grey-levelvalue of the instantaneous pixel exceeds the threshold value T_(i-1)(the reflectivity of the wall exceeds that of the blood). If it happensthat an edge of the image interrupts the mark of the vessel, theestimated end of the vessel in this direction does not constitute a truewall, but still the pixels between two angularly adjacent rays willconstitute a fraction of the local mark of the vessel. The latter stepprovides sampling of segments forming the limits of the vessel in theimage and, by joining the extremities of each pair of successive rays, arepresentation of the vessel is obtained which is composed of triangularsectors 12bis as shown in FIG. 4.

These sectors are characterized by the coordinates of their externalapexes P1_(j), P1_(j+1), the third point being the common origin Pi. Foreach sector six analytical formulas are applied to the coordinates ofits two external apexes so as to supply the local information necessaryto enable ultimately a linear regression to be carried out on the entirelocal mark of the vessel.

The value of the Doppler angle DA (FIG. 5) is thus deduced directly fromthe slope a of the regression line 13, the direction of the ultrasonictransmission axis (or excitation axis) which is known (vector 10), andfinally a correlation coefficient r which can be calculated in paralleland supplies an indication as regards the relevance of the resultobtained.

The execution of the entire process is summarized in the organigram ofFIG. 6. During the first step, at 15, the coordinates of the initialpoint (or pixel) Pi are stored in the memory 4 (see FIG. 1). At 16, thefirst ray tracing is performed, starting from the initial point. At 17,the threshold T_(i-1) is determined on the basis of the histogram ofselected points. At 18, a second ray tracing, associated with thethreshold T_(i-1), enables formation of the first local mark which isconstituted by assembled triangular sectors 12. At 19, the calculationsand accumulations for the sectors are performed. At 20, the resultsobtained during the preceding step are used for the linear regressionstep (calculation of a and of r), and finally the result, being theangle DA, is obtained during step 21.

The following is an example of the execution of the described process.It will be noted first of all that the process considered requires thesame sampling step on the horizontal and vertical axes; therefore, itmust be applied to images which have been subjected to a digital scanconversion (DSC) operation. As DSC is a prior necessity for treatment,the images may be acquired from any type of transducer network.Actually, the image used as a basis is the reflectivity image in greylevels which is customarily displayed on the screen of the system (ofthe workstation).

RAY TRACING AND FORMATION OF THE HISTOGRAM

First of all, the operator performs a pointing operation in the image bydesignation of a location in the vessel. The ray tracing operation thusconsists in the recursive addition (by simple programming which will beevident to those skilled in the art) of pairs of shifts in x and in y(each pair characterizing a ray) to the coordinates of this locationuntil the edge of the image is reached. For each intermediate position apixel of the image is addressed along a linear ray (see FIG. 3).

Evidently, the larger the number of rays, the better the sampling of thearea of the vessel will be. In practice it has been found that takingless than 32 rays could occasionally be insufficient, whereas choosing32 rays has always given good results with the experimentally obtaineddata. Therefore, it is suggested to use preferably 32 rays.

If the number of grey levels in the image is denoted as NGL, thisprocessing step supplies a set of NGL integers, each value of whichindexed by q (q=0, 1, . . . , NGL), has been incremented by 1 for eachpixel of a value indexed by q along the course of the rays.

ISODATA COALESCENCE ALGORITHM

This algorithm aims to determine threshold values so as to classify thepixels of the image segmented in nc groups having a physicalsignificance in respect of the contents of the image.

In this case 3 (nc=3) classes (VESSEL, INTERMEDIATE, WALL) are chosenwhose behaviour (the reflectivity) is well known, and the thresholdvalue T₂ between the classes INTERMEDIATE and WALL is needed. Thisgrouping algorithm performs an overall segmentation of the image bydetermination of threshold values deduced from the histogram of the greylevels of the image which, therefore, is not dependent (at least not inthe first order) on the location of the object to be segmented. Thisalgorithm necessitates only an initial choice of the value of the meangrey level for each class, and iterative updating of these values isterminated by the supply of threshold values between classes.

First Initialization Step:

Let h(j) be the probability density of the grey-level value j in theoriginal image.

Let min,max! be the smallest grey-level value interval encompassing thevalues of h(j) which are not zero.

Let m_(i), i ε{1, . . . , nc} be the mean initial value for all classes.This initial estimation may be carried out by dividing the axis of thegrey levels into nc equidistant intervals and by calculating thearithmetical mean value over each interval as if the probability densitywas uniform between min and max.

Second Step:

The (nc-1) threshold values T_(i) are estimated by utilizing thefollowing relation (where x! represents the integer pan of x):

    T.sub.i = (m.sub.i.sbsb.-- m.sub.i+1)/2!,i ε{1 . . . , nc-1}

All pixels whose grey-level values are situated in the interval:

    A.sub.i = T.sub.i-1, T.sub.i !,i ε{1, . . . , nc}

are then assigned to the class i (where T_(o) =min-1 and T_(nc) =max).

Third Step:

The mean values of the classes are updated by utilizing the followingrelation (where x! represents the integer part of x): ##EQU7##

Fourth and Last (Iteration) Step:

If at least one of the m_(i) values has been modified during the thirdstep, it is necessary to return to the second step, thus performing aloop until the algorithm converges; if not, the final threshold valuesare obtained. This processing step produces an integer threshold valueT_(i-1) which separates the blood from the arterial walls across theimage.

WALL DETECTION AND REPRESENTATION OF THE VESSEL IN SECTORS

Once the desired threshold value T_(i-1) has been obtained, the same raytracing method is carried out as described above in order to detect thearea in which the wall of the vessel commences, along each ray scannedsimply, by comparing the values of pixels along a ray with the thresholdvalue. The coordinates of the last pixel marking the presence of thewall, along each ray, is stored in the memory (the memory 4 in FIG. 1).

Thus, a sampling of the limits of the wall of the vessel in the image isobtained and, by joining each pair of ends of rays which succeed oneanother angularly, a representation of the vessel is obtained which iscomposed of triangular radial sectors (see FIG. 4).

Let nr be the number of rays. This processing step produces two sets ofnr integer numbers whose values, indexed by s, correspond to thecoordinates of the last point of the ray indexed by s, so Pl_(s).

CHARACTERIZATION OF THE REPRESENTATION IN SECTORS

The actual object of the whole operation is to supply the following twoparameters:

The slope a of the regression line deduced by the conventionalleast-squares method, supplying the tangent of the Doppler angle DA(DA=Arc tg a): ##EQU8## where {x(n),y(n)}, n ε {1, . . . , N}, are thecoordinates of the pixels comprised in the representation of the vesselin triangular sectors, i.e. pixels relating to the local mark.

The correlation coefficient r whose absolute value will be nearer to 1as the representation in sectors is nearer to a collection of pointspresenting a regular alignment along a given axis (0≦|r|≦1): ##EQU9##From a practical point of view, the determination of these two values (aand r) requires measurement and accumulation, over the whole local mark,of six different quantities (six terms) deriving a contribution from thecoordinates of each pixel contained therein: ##EQU10## The justificationof the decomposition of the local mark of the vessel into triangularsectors is found here: the values of the above parameters,characterizing all points contained in a given sector, can indeed bededuced from the coordinates of the two external apexes of this sectorexclusively. This signifies that the characterization of a sector isrealised by application of 6 analytical formulas to 4 integer values asexplained hereinafter.

Let xm,y(xm)! and xM,y(xM)! be the coordinates of the two externalapexes (Pl_(i),Pl_(i+1)). where xM≧xm, the sector considered beingsituated by symmetry, in the first trigonometric quadrant from the topand the right (by multiplication of the negative coordinates by -1 untilx>0 and y>0). 6 types of different possible sectors are obtained, beingin this configuration restricted to the first quadrant, in dependence onthe values and the relative sequencing of the coordinates of theexternal apexes as shown in FIG. 7 (see notably line 76 of FIG. 7).

In this FIG. 7 classes of sectors can be distinguished which representall possible shapes of sectors, said classes being shown as 7 lineswhich are referenced as 71, 72, 73, 74, 75, 76, 77 from the top down,the lines 71, 72 and 73 each representing two sector shapes which bothrelate to the same respective class. These classes are characterized asfollows:

71: y(xm)=y(xM)

72: xm=xM

73: xm=0

74: y(xM)=0

75: y(xm)=0

76: y(xm)>y(xM) x xm/xM

77: y(m)<y(xM) x xm)xM.

The classes of the lines 76 and 77, being characterized by inequalitiesin the opposite sense in relation to the same expression, may beconsidered as belonging to the same class from a mathematical point ofview.

The above-mentioned analytical formulas used for the calculations are:##EQU11## where: f(x,y)=1 for deducing the pixels and obtaining N

f(x,y)=x in order to obtain sx

f(x,y)=y for obtaining sy

f(x,y)=xy for obtaining sxy

f(x,y)=xx for obtaining sx2

f(x,y)=yy for obtaining sy2.

The integration limits of course, are dependent on the shape of thesector as indicated in Annex 1 (relating to the classes of FIG. 7). butthe ultimate expressions to be calculated so as to perform the summingfor each type of sector are very simple and do not require a substantialcalculation effort as indicated in annex 2 (relating to the classes ofFIG. 7).

When each sector of the representation of the vessel has been treated,the actual regression can be performed.

LINEAR REGRESSION AND ESTIMATION OF THE VALIDITY OF THE RESULT

The Doppler angle DA and the reliability criterion r are obtaineddirectly from 6 values supplied at the end of the foregoing step, aftertreatment of all sectors: ##EQU12## On the other hand: ##EQU13## Thisvalue r is a control parameter which may be used to form a reliabilitycriterion for the measurement: if its value is larger than, for example0.5, the result may be considered to be reliable; if not, there isproblem which prevents the system from supplying a valid result (notablya local mark having a circular shape due to poor positioning of thesectional plane) and the operator will have to make another attemptafter correction.

The basic method described above can be refined by applying it tosuccessive areas within the vessel by automatic centring.

Once a first representation in sectors of the vessel has been realised,it is easy to calculate the coordinates of its centre of gravity G (forexample, by taking the mean value of the abscissae and the ordinates ofthe pixels of the first local mark and, while choosing this point G asthe new initial point, by repeating the entire process described aboveso as to determine a second local mark which has been optimized and alsoan optimized Doppler angle). This enhancement of the method enables alarger margin in respect of the choice of the initial point so that theuser can benefit from automatic centring within the vessel (for example,it is thus possible to indicate, using the mouse of the workstation, alocation close to a wall without being penalized), and the resultsobtained will be more reliable (standard deviation of less than 1°).

It is to be noted that this first repeat of the entire process may bereiterated, for example a predetermined number of times which may be ashigh as 20, or until convergence on a particular point is obtained. Thelatter point corresponds to the centre of gravity of the entire trace ofthe vessel in the echographic image and not only to that of the firstlocal mark deduced from the choice of the first initial point by theuser. Thus, the method no longer concerns a local measurement, butoffers the best environment along the entire trace of the vessel inorder to extract the value of the Doppler angle therefrom.

    ______________________________________                                        ANNEX 1: Integration limits for each triangular sector configuration          ______________________________________                                        line 1:      a.sub.0 = y(xm)/xM                                               y(xm) = y(xM)                                                                              a.sub.1 = y(xm)/xm                                                            a.sub.2 = y(xm)/xM                                                            b.sub.2 = a.sub.3 = 0                                                         b.sub.3 = y(xm)                                                  line 2:      a.sub.0 = MIN y(xm), y(xM)!/xm                                   xM = xm      a.sub.1 = MAX y(xm), y(xM)!/xm                                                a.sub.2 = b.sub.2 = a.sub.3 = b.sub.3 = 0                        line 3:      a.sub.0 = a.sub.1 = b.sub.2 = 0                                  xm = 0       a.sub.2 = y(xM)/xM                                                            a.sub.3 =  y(xM) - y(xm)!/xM                                                  b.sub.3 = y(xm)                                                  line 4:      a.sub.0 = a.sub.2 = b.sub.2 = 0                                  y(xM) = 0    a.sub.1 = y(xm)/xm                                                            a.sub.3 = -y(xm)/xM - xm                                                      b.sub.3 = y(xm) - a.sub.3 × xm                             line 5:      a.sub.0 = b.sub.3 = 0                                            y(xm) = 0    a.sub.1 = y(xM)/xM                                                            a.sub.2 = y(xM)/(xM - xm)                                                     b.sub.2 = -a.sub.2 × xm                                                 a.sub.3 = y(xM)/xM                                               lines 6 and 7:                                                                             a.sub.0 = y(xm)/xm or y(xM)/xM                                   y(xm)        a.sub.1 = y(xM)/xM or y(xm)/xm                                   < or >       a.sub.2 =  y(sM) - y(xm)!/(xM - xm) or y(xM)/xM                  y(xM) × xm/xM                                                                        b.sub.2 = y(xm) - a.sub.2 × xm or 0                                     a.sub.3 = y(xM)/xM or  y(xM) - y(xm)!/(sM - xm)                               b.sub.3 = 0 or y(xm) - a.sub.3 × xm                        ______________________________________                                    

    __________________________________________________________________________    ANNEX 2: Expression to be calculated (in this table s += x signifies: add     x to s)                                                                       __________________________________________________________________________    line 1:  sn += y(xm) × (xM - xm)/2                                      y(xm) = y(xM)                                                                          sx += y(xm) × (xM) - xm) × (xM + xm)/6                            sy += y(xm).sup.2 × (xM) - xm)/3                                        sxy += y(xm).sup.2 × (xM - xm) × (xM + xm)/8                      sx2 += y(xm) × (xM.sup.3 - xm.sup.3)/12                                 sy2 += y(xm).sup.3 × (xM - xm)/4                               line 2:  YMAX = MAX y(xm), y(xM)!                                             xM = xm  YMIN = MIN y(xm), y(xM)!                                                      sn += xm × (YMAX - YMIN)/2                                              sx += xm.sup.2 × (YMAX - YMIN)/3                                        sy += xm × (YMAX - YMIN) × (YMAX + YMIN)/6                        sxy += xm.sup.2 × (YMAX - YMIN) × (YMAX + YMIN)/8                 sx2 += xm.sup.3 × (YMAX - YMIN)/4                                       sy2 += xm × (YMAX.sup.3 - YMIN.sup.3)/12                       line 3:  sn += xM × y(xm)/2                                             xm = 0   sx += xM.sup.2 × y(xm)/6                                                sy += xM × y(xm) ×  y(xM) + y(xm)!/6                              sxy += xM.sup.2 × y(xm) ×  2 × y(xM) +                      y(xm)!/24                                                                     sx2 += xM.sup.3 × y(xm)/12                                              sy2 += xM × y(xm) ×  y(xM).sup.2 + y(xM) ×                  y(xm) + y(xm).sup.2 !/12                                             line 4:  sn += xM × y(xm)/2                                             y(xM) = 0                                                                              sx += xM × y(xm) × (xM + xm)/6                                    sy += xM × y(xm).sup.2 /6                                               sxy += xM × y(xm).sup.2 × (xM + 2 × xm)/24                  sx2 += xM × y(xm) × (xM.sup.2 + xM × xm +                   xm.sup.2)/12                                                                  sy2 += xM × y(xm).sup.3 /12                                    line 5:  sn += xm × y(xM)/2                                             y(xm) = 0                                                                              sx += xm × y(xM) × (xM + xm)/6                                    sy += xm × y(xM).sup.2 /6                                               sxy += xm × y(xM).sup.2 × (2 × xM + xm)/24                  sx2 += xm × y(xM) × (xM.sup.2 + xM × xm +                   xm.sup.2)/12                                                                  sy2 += xm × y(xM).sup.3 /12                                    lines 6 and 7:                                                                         sn +=  xm × y(xM) - xM × y(xm!/2                         y(xm) <  sx += (xM + xm) ×  xm × y(xM) - xM × y(xm)!/6      y(xM) × xm/xM                                                                    sy +=  y(xM) + y(xm)!                                                           ×  xm × y(xM) - xM × y(xm)!/6                    if       sxy +=  xm × y(xM) - xM × y(xm)!                         y(xm) >    ×  2 × xM × y(xM) + xm × y(xM)             y(xM) × xm/xM                                                                        + xM × y(xm) + 2 × xm × y(xm)!/24              multiply each                                                                          sx2 +=  xm × y(xM) - xM × y(xm)!                         expression by -1.                                                                        × (xM.sup.2 + xM × xm + xm.sup.2)/12                            sy2 +=  xm × y(xM) - xM × y(xm)!                                    ×  y(xM).sup.2 + y(xM) × xm) + y(xm).sup.2             __________________________________________________________________________             !/12                                                             

I claim:
 1. A method for automatic determination of Doppler angle DA,enclosed by the direction of an echographic excitation and the axis of arelevant blood vessel in an ultrasonic echographic grey-level image, onthe basis of prior designation of an initial point situated in thevicinity of the axis of the relevant vessel, comprising the followingsteps:a) first isotropic tracing of rays from said initial point whichtraverse the entire echographic image and forming a histogram of greylevels of selected points which are regularly distributed along saidrays, b) applying an image processing algorithm to said histogram inorder to classify the grey levels of said selected points into at leasttwo classes CL₁, CL₂, . . . , two adjacent classes being separated bythresholds T₁, T₂, . . . , expressed in grey levels, one of the classesCL_(i), bounded by a lower threshold T_(i-1), being representative ofwalls of blood vessels in said image, c) second tracing of rays fromsaid initial point, during which the grey level of each point of eachray is compared with the threshold T_(i-1) and each ray is restricted tothe first end point Pl encountered whose grey level is equal to orhigher than T_(i-1), and producing a representation of said relevantblood vessel in the form of a first local mark which is composed oftriangular sectors which have said starting point within said firstlocal mark in common, d) determining a slope (a) of the regression lineof said first local mark by applying a linear regression method to the Npixels, having the coordinates x(n) and y(n) of said first local mark,as: ##EQU14## e) and calculating a Doppler angle as:

    DA=Arc tg (a).


2. A method as claimed in claim 1, further comprising a step ofvalidating the value of the Doppler angle DA only for a correlationcoefficient r of the pixels of said first local mark which exceeds athreshold Rm of predetermined value between 0 and 1, where: ##EQU15## 3.A method as claimed in claim 1, further comprising a step of validatingthe value of the Doppler angle DA only if the symmetry of said localmark with respect to its regression line (or its inertial axis) isconsidered to be sufficient, said symmetry being deduced from thecalculation of a centered moment of order 3, being sk(f): ##EQU16## andf is the function:

    y(n)=f(x(n)), where n is within  1, . . . , N!,

the centered moment of order 3 being calculated for two curves which arerepresentative of the function f: sk1 for the curve formed by theextremities of rays of said local mark which are situated above theregression line (or the inertial axis) taken as the axis of theabscissae (y>o), sk2 for the curve formed by the extremities of rayssituated below aid line (y<o), the values of sk1 and sk2 being comparedafterwards.
 4. A method as claimed in claim 1, characterized in that theterms for calculating the slope (a) of the regression line and of thecorrelation coefficient r are calculated for each triangular sectorformed by two rays whose ends have the coordinates xm, y(xm) and xM,y(xM), on the basis of the following analytical formula: ##EQU17## inwhich the function f(x,y) is given the following successive values: 1;x; xy; xx; yy, and the coefficients a₀, a₁, a₂, b₂, a₃, b₃ are givensets of values which are dependent only on the coordinates of said endpoints of two rays, with 6 sets of possible values associated with 6respective different classes of different triangular sectors.
 5. Amethod as claimed in claim 1, further comprising a step for positioningwhich comprises calculating a center of gravity G of said first localmark, said center of gravity G then being chosen as a new initial pointwhich has been optimized for subsequent repetition of the method forautomatic determination of Doppler angle DA.
 6. A method for measurementas claimed in claim 5, wherein said step for positioning and thesubsequent repetition of the method for automatic determination ofDoppler angle DA are repeated a predetermined number of times which isbetween 1 and
 20. 7. An arrangement for automatic determination ofDoppler angle DA, enclosed by the direction of an echographic excitationand the axis of a blood vessel in an ultrasonic echographic grey-levelimage, on the basis of prior designation of a initial pixel situated inthe vicinity of the axis of the relevant vessel, which arrangementcomprises:a memory for storing the grey-level values of the echographicimage in the form of a matrix of pixels, a screen for displaying theechographic image in the form of a two-dimensional matrix of pixels,pointer means for pointing out said initial pixel in the displayedimage, and calculation means which can access said memory and saidinitial point,the calculation means being arranged to perform: firstisotropic tracing of rays from said initial pixel which traverse theentire echographic image and formation of a histogram of grey levels ofpixels which are selected so as to be regularly distributed along saidrays, application of an image processing algorithm to said histogram inorder to classify the grey levels of said selected pixels in at leasttwo classes CL₁, CL₂, . . . , two adjacent classes being separated bythresholds T₁, T₂, . . . , expressed in grey levels, one of the classesCL_(i), bounded by the lower threshold T_(i-1), being representative ofwalls of blood vessels in said image, second tracing of rays from saidinitial pixel, during which the grey level of each pixel of each ray iscompared with the threshold T_(i-1) and each ray is restricted to thefirst end pixel Pl encountered whose grey level is equal to or higherthan T_(i-1), and production of a representation of said blood vessel inthe form of a first local mark which is composed of triangular sectorswhich have said starting pixel within said first local mark in common,determination of a slope (a) of the regression line of said first localmark by application of a linear regression method to the N pixels,having the coordinates x(n) and y(n) of said first local mark, as:##EQU18## and calculation of a Doppler angle as:

    DA=Arc tg (a).


8. A method for automatic measurement as claimed in claim 1, comprisinga step of automatically displaying the value of said Doppler angle DA inthe field of a screen which contains the echographic image in the formof a two-dimensional matrix of pixels in grey levels and/or color, andof drawing a line segment representing the axis of the relevant vessel.9. A method for automatic measurement as claimed in claim 1, whereinspeed information is represented in the form of an image, or color speedmap, with at least two colors, each of which is a function of the signof the speed and a color intensity proportional to the absolute value ofthe speed, a speed color scale being represented in a field of thescreen provided for this purpose, comprising a step of automaticallycorrecting said speed color scale, so as to ensure that it contains adynamic of the speeds observed along the axis of said relevant bloodvessel, by automatically displaying the maximum positive and negativevalues of these speeds at the two extremities of said scale.
 10. Amethod for measurement as claimed in claim 1, wherein evolution of aspeed distribution in the course of time in a segment of the relevantvessel is represented in a mode M, on a screen provided for thispurpose, in the form of a histogram showing a dynamic of speeds as afunction of time, referred to as a Doppler spectrum, comprising a stepof automatically correcting the scale of the speeds represented in saidDoppler spectrum in order to ensure that it contains the dynamic of thespeeds observed along the axis of said blood vessel.
 11. A method formeasurement as claimed in claim 1, further comprising a step ofvalidating the value of the Doppler angle DA only for a correlationcoefficient r of the pixels of said first local mark which exceeds athreshold Rm of predetermined value between 0 and 1, where: ##EQU19##